home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libconv / TRY / ornl_dfir2d.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  4.1 KB  |  134 lines
  1.       subroutine ornl_dfir2d(    f,incf,ldf,ifx0,n_fx,ify0,n_fy,
  2.      $                g,incg,ldg,igx0,n_gx,igy0,n_gy,
  3.      $                h,inch,ldh,ihx0,n_hx,ihy0,n_hy,
  4.      $                alpha, beta)
  5.       double precision f(ldf,*)
  6.       double precision g(ldg,*)
  7.       double precision h(ldh,*)
  8.       double precision alpha, beta
  9.       integer incf, ldf, ifx0, n_fx, ify0, n_fy
  10.       integer incg, ldg, igx0, n_gx, igy0, n_gy
  11.       integer inch, ldh, ihx0, n_hx, ihy0, n_hy
  12. c-----------------------------------------------------------------------------
  13. c
  14. c   dfir2d  performs a 2D convolution in the time domain :
  15. c    h(i,j) = beta * h(i,j) + alpha * Sum.Sum[ f(k,l) * g(i-k,j-l) ]
  16. c
  17. c-----------------------------------------------------------------------------
  18. c
  19. c   PARAMETERS:
  20. c
  21. c    f:    Vector containing sequence "f"
  22. c    incf:    Increment between two successive lines   of "f"
  23. c    ldf:    Increment between two successive columns of "f"
  24. c    ifx0:    Index of the first element of each column of "f"
  25. c    n_fx:    Number of elements (lines) of each column of "f"
  26. c    ify0:    Index of the first column of "f"
  27. c    n_fy:    Index of the  last column of "f"
  28. c    
  29. c    g:    Vector containing sequence "g"
  30. c    incg:    Increment between two successive lines   of "g"
  31. c    ldg:    Increment between two successive columns of "g"
  32. c    igx0:    Index of the first element of each column of "g"
  33. c    n_gx:    Number of elements (lines) of each column of "g"
  34. c    igy0:    Index of the first column of "g"
  35. c    n_gy:    Index of the  last column of "g"
  36. c    
  37. c    h:    Vector containing sequence "h"
  38. c    inch:    Increment between two successive lines   of "h"
  39. c    ldh:    Increment between two successive columns of "h"
  40. c    ihx0:    Index of the first element of each column of "h"
  41. c    n_hx:    Number of elements (lines) of each column of "h"
  42. c    ihy0:    Index of the first column of "h"
  43. c    n_hy:    Index of the  last column of "h"
  44. c
  45. c    alpha: Scaling factor for the Convolution
  46. c    beta:  Scaling Factor for the Output on Entry
  47. c
  48. c-----------------------------------------------------------------------------
  49. c
  50. c   PLEASE NOTE:
  51. c    Please note that the array pointers must all point to the first 
  52. c    element of the array "(ifx0,ify0)", "(igx0,igy0)" and "(ihx0,ihy0)"
  53. c     If "f" for example is defined as
  54. c        dimension f(-25:45,10:21)
  55. c    Then "dfir2d" must be called with the following parameters
  56. c        call dfir2d( f(-25,10),(45-(-25)+1),-25,45,10,21 ... )
  57. c
  58. c-----------------------------------------------------------------------------
  59. c
  60. c   This Fortran Subroutine written by
  61. c    Jean-Pierre Panziera
  62. c    Silicon Graphics Inc
  63. c    September 27, 1991
  64. c
  65. c-----------------------------------------------------------------------------
  66.  
  67.       double precision zero, one
  68.       parameter (zero = 0.0, one = 1.0)
  69.       integer iy0, iy1, ify1, igy1, ihy1
  70.       integer k1, k2, k, j, ix, iy, jf, jg, jh
  71.  
  72. c-----------------------------------------------------------------------------
  73. c
  74. c   Compute Y Boundaries
  75. c
  76. c-----------------------------------------------------------------------------
  77.     ify1 = ify0 + n_fy - 1
  78.     igy1 = igy0 + n_gy - 1
  79.     ihy1 = ihy0 + n_hy - 1
  80.  
  81.     iy0 = ihy0
  82.     if( iy0 .lt. (ify0+igy0))    iy0 = ify0+igy0
  83.     iy1 = ihy1
  84.     if( iy1 .gt. (ify1+igy1))    iy1 = ify1+igy1
  85.  
  86. c-----------------------------------------------------------------------------
  87. c
  88. c   Scale the output
  89. c
  90. c-----------------------------------------------------------------------------
  91.       if( beta .eq. zero) then
  92.     do j = 1, n_hy
  93.         ix = 1
  94.         do k = 1, n_hx
  95.         h(ix,j) = zero
  96.         ix = ix + inch
  97.         end do
  98.     end do
  99.       else if( beta .ne. 1.0) then
  100.     do j = 1, n_hy
  101.         ix = 1
  102.         do k = 1, n_hx
  103.         h(ix,j) = beta * h(ix,j)
  104.         ix = ix + inch
  105.         end do
  106.     end do
  107.       endif
  108.  
  109.       if( alpha .eq. zero ) return
  110.  
  111. c-----------------------------------------------------------------------------
  112. c
  113. c   Call the 1D convolution 
  114. c
  115. c-----------------------------------------------------------------------------
  116.     do j = iy0, iy1
  117.         k1 = max( igy0, j-ify1)
  118.         k2 = min( j-ify0, igy1)
  119.         jh = j-ihy0+1
  120.         do k = k1, k2
  121.         jf = j-k-ify0+1
  122.         jg = k-igy0+1
  123.         call dfir1d(    f(1,jf), incf, ifx0, n_fx,
  124.      $                g(1,jg), incg, igx0, n_gx,
  125.      $                h(1,jh), inch, ihx0, n_hx,
  126.      $                alpha, one)
  127.         end do
  128.     end do
  129. c-----------------------------------------------------------------------------
  130.       return
  131.       end
  132.  
  133.  
  134.